MODELLI LINEARI¶
I modelli lineari utilizzati saranno modelli ARIMA e modelli UCM.
Nella maggior parte delle aziende manufatturiere, l'analisi delle serie storiche guida le attività fondamentali di pianificazione aziendale, approvvigionamento e produzione. Eventuali errori nelle previsioni si propagheranno lungo tutta la catena di fornitura o in qualsiasi contesto aziendale. Quindi è importante ottenere previsioni accurate per risparmiare sui costi ed è fondamentale per il successo.
Non solo nella produzione, le tecniche ed i concetti alla base della previsione delle serie temporali sono applicabili in qualsiasi azienda.
Ora la previsione di una serie temporale può essere sostanzialmente suddivisa in due tipi.
Se si utilizzano solo i valori precedenti della serie storica per prevederne i valori futuri, si parla di previsione di serie storica univariata.
E se si utilizzano predittori diversi dalle serie (note anche come variabili esterne) per fare previsioni, si parla di previsione multivariata delle serie temporali.
ARIMA, abbreviazione di "AutoRegressive Integrated Moving Average", è un algoritmo di previsione basato sull'idea che le informazioni nei valori passati delle serie temporali possono essere utilizzate da sole per prevedere i valori futuri.
Caricamento delle librerie, path e dataset
# import necessary libraries
import os
import pandas as pd
import re
import numpy as np
import plotly.graph_objects as go
import plotly.express as px
import plotly.io as pio
# Themes
# "seaborn", 'ggplot2', 'seaborn', 'simple_white', 'plotly', 'plotly_white', 'plotly_dark', 'presentation', 'xgridoff', 'ygridoff', 'gridon', 'none'
pio.templates.default = "presentation"
from plotly.subplots import make_subplots
# select data folder_dati
glob_path = os.path.normpath(os.getcwd() + os.sep + os.pardir)
folder_dati = glob_path+R"//dati"
folder_risultati = glob_path+R"//risultati"
files = os.listdir(folder_risultati)
files[0:]
# Load export.csv, original file
df = pd.read_csv(folder_risultati+"/export.csv", sep = ",")
# substitute null values with 0
df['Minuti ciclo'] = df['Minuti ciclo'].fillna(0)
df['Minuti NO ciclo'] = df['Minuti NO ciclo'].fillna(0)
df['Efficienza'] = df['Efficienza'].fillna(0)
# convert Date
df['Date'] = pd.to_datetime(df["Date"], format='%Y-%m-%d')
# convert minuti ciclo, no ciclo, efficienza to double
df['Minuti ciclo'] = df['Minuti ciclo'].replace(',', '.', regex=True).astype(float)
df['Minuti NO ciclo'] = df['Minuti NO ciclo'].replace(',', '.', regex=True).astype(float)
df['Efficienza'] = df['Efficienza'].replace(',', '.', regex=True).astype(float)
df_grouped_date = df.groupby(['Date'])['Minuti ciclo', 'Minuti NO ciclo'].sum().reset_index()
df_grouped_date['Efficienza'] = df_grouped_date['Minuti ciclo'] / (df_grouped_date['Minuti ciclo'] +
df_grouped_date['Minuti NO ciclo'])
###################### CLIENTE
cliente = 'CAMA'
df_cliente = df
df_cliente = df_cliente[df_cliente['CLIENTE'] == cliente]
df_grouped_date_cliente = df_cliente.groupby(['Date'])['Minuti ciclo', 'Minuti NO ciclo'].sum().reset_index()
# create field date
df_days = pd.DataFrame({'Date': df_grouped_date.Date})
merged = pd.merge(df_days, df_grouped_date_cliente, on=['Date'], how='left')
df_grouped_date_cliente = merged.fillna(0)
# days per week
eff_week = df_grouped_date_cliente
eff_week['Giorno'] = eff_week['Date'].dt.strftime("%A")
df_grouped_date_cliente = eff_week
# delete sundays
df_grouped_date_cliente = df_grouped_date_cliente[df_grouped_date_cliente['Giorno'] != "Sunday"]
# i valori a 0 sono un po... sostituisco con media
df_grouped_date_cliente['Minuti ciclo'] = df_grouped_date_cliente['Minuti ciclo'].mask(df_grouped_date_cliente['Minuti ciclo']==0).fillna(df_grouped_date_cliente['Minuti ciclo'].mean())
df_grouped_date_cliente.to_csv(folder_risultati + '//df_grouped_date_cliente.csv', decimal = ',', index = False)
df_grouped_date_cliente
df_grouped_risorsa = df.groupby(['Date', 'Risorsa'])['Minuti ciclo', 'Minuti NO ciclo'].sum().reset_index()
df_grouped_risorsa['Efficienza'] = df_grouped_risorsa['Minuti ciclo'] / (df_grouped_risorsa['Minuti ciclo'] +
df_grouped_risorsa['Minuti NO ciclo'])
df_grouped_date.to_csv(folder_risultati + '\\df_grouped_date.csv', decimal = ',', index = False)
df_grouped_risorsa.to_csv(folder_risultati + '\\df_grouped_date_risorsa.csv', decimal = ',', index = False)
grouped_vtc530 = df_grouped_risorsa[df_grouped_risorsa['Risorsa'] == "MAZAK-VTC530"]
grouped_vtc800 = df_grouped_risorsa[df_grouped_risorsa['Risorsa'] == "MAZAKVTC800"]
grouped_j600= df_grouped_risorsa[df_grouped_risorsa['Risorsa'] == "MAZAKJ600"]
grouped_2014 = df_grouped_risorsa[df_grouped_risorsa['Risorsa'] == "MAZAK2014"]
grouped_2015 = df_grouped_risorsa[df_grouped_risorsa['Risorsa'] == "MAZAK2015"]
# Now divide also by hour
df['Hour'] = df['KEY'].str.rpartition('_')[2]
df.head()
# new format of date
df_hour = df
df_hour['Date'] = pd.to_datetime(df['KEY'].str.split('_').str[0] + "/" +
df['KEY'].str.split('_').str[1] + "/" +
df['KEY'].str.split('_').str[2] + " " +
df['Hour'] + ":00:00", format='%Y/%m/%dT%H:%M:%S')
df_hour.head()
df_grouped_hour = df_hour.groupby(['Date'])['Minuti ciclo', 'Minuti NO ciclo'].sum().reset_index()
df_grouped_hour['Efficienza'] = df_grouped_hour['Minuti ciclo'] / (df_grouped_hour['Minuti ciclo'] +
df_grouped_hour['Minuti NO ciclo'])
# Efficiency per week
eff_week = df_grouped_date
eff_week['Giorno'] = eff_week['Date'].dt.strftime("%A")
df_grouped_date = eff_week
fig = px.box(eff_week, x="Giorno", y="Efficienza", color="Giorno")
fig.update_traces(quartilemethod="exclusive") # or "inclusive", or "linear" by default
fig.show()
Dal grafico soprastante si evince che mediamente l'efficienza delle macchine cresce al crescere dei giorni nella settimana.
# delete sundays
df_grouped_date = df_grouped_date[df_grouped_date['Giorno'] != "Sunday"]
df_grouped_date.to_csv(folder_risultati + '\\df_grouped_date_nosunday.csv', decimal = ',', index = False)
# plot general efficiency
fig = px.line(df_grouped_date_cliente, x = "Date", y = "Minuti ciclo",
title = "minuti ciclo cliente / gg")
fig.show()
Grafici e visualizzazioni¶
I grafici saranno della serie di dati raggruppata per giorno e per ora
# plot general efficiency
fig = px.line(df_grouped_date, x = "Date", y = "Efficienza",
title = "Efficienza / giorno Totale delle macchine")
fig.show()
Dal primo grafico si può notare il trend crescente che ha portato l'efficienza media delle macchine dal 30% al 45% nel giro di un anno e mezzo!
Per quanto riguarda la stagionalità , è presente nel senso che la domenica le macchine non lavorano, mentre al sabato lavorano soltanto 5/6 ore: dunque si ha una stagionalità settimanale.
# plot general efficiency per hour
fig = px.line(df_grouped_hour, x = "Date", y = "Efficienza",
title = "Efficienza / giorno Totale delle macchine")
fig.show()
Anche dal grafico orario si intravede un trend crescente da inizio 2019 fino a metà 2020.
In questo caso la stagionalità è rappresentata prima dai periodi 6-22 di lavoro della macchina, poi dalla stagionalità settimanale prima citata.
In questo grafico vengono evidenziate meglio i due periodi di chiusura aziendale, coincidenti con la chiusura natalizia (21 dicembre - 2 gennaio compresi) e la chiusura causa pandemia da COVID-19 (19 - 31 marzo compresi)
# plot general efficiency per day and macchine
fig = px.line(df_grouped_risorsa, x = "Date", y = "Efficienza", color = "Risorsa",
title = "Efficienza / giorno Totale delle macchine")
fig.show()
# plot Minuti ciclo general
fig = px.line(df_grouped_date, x = "Date", y = "Minuti ciclo",
title = "Minuti ciclo / giorno")
fig.show()
# plot Minuti NO ciclo general
fig = px.line(df_grouped_date, x = "Date", y = "Minuti NO ciclo",
title = "Minuti NO ciclo / giorno")
fig.show()
# save
df_grouped_hour.to_csv(folder_risultati + '\\df_grouped_hour.csv', decimal = ',', index = False)
df_hour.to_csv(folder_risultati + '\\df_hour.csv', decimal = ',', index = False)
Viene effettuata una decomposizione delle componenti della serie storica
from statsmodels.tsa.seasonal import seasonal_decompose
# crea la serie storica
ts = df_grouped_date[['Date', 'Minuti ciclo']]
ts.Date = pd.to_datetime(ts.Date)
ts.set_index('Date', inplace=True, drop=True)
ts['Minuti ciclo'] = ts['Minuti ciclo'].fillna(0)
decompose_efficiency = seasonal_decompose(ts['Minuti ciclo'], model='additive', freq=7)
print(ts.info())
# Plot
decompose_efficiency.plot().show()
Le componenti decomposte mostrano un trend comunque crescente, con una stagionalità settimanale.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime
from datetime import timedelta
from pandas.plotting import register_matplotlib_converters
from statsmodels.tsa.stattools import acf, pacf
from statsmodels.tsa.statespace.sarimax import SARIMAX
register_matplotlib_converters()
from time import time
# esporto solo i dati che mi interessano
ts_f = ts[:491]
# i valori a 0 sono 57, contando nnn gg di chiusura aziendale + le varie festività nel giro di due anni ci sta! Li sostituisco con la media!
ts_f['Minuti ciclo'] = ts_f['Minuti ciclo'].mask(ts_f['Minuti ciclo']==0).fillna(ts_f['Minuti ciclo'].mean())
# train e test.....
train_end = datetime(2019,12,31)
test_end = datetime(2020,7,23)
train = ts_f[:train_end]
test = ts_f[train_end + timedelta(days=1):test_end]
MODELLO ARIMA¶
ARIMA, abbreviazione di "Auto Regressive Integrated Moving Average" è in realtà una classe di modelli che spiega una data serie temporale in base ai propri valori passati, ovvero i propri ritardi e gli errori di previsione ritardati, in modo che sia possibile utilizzare l'equazione per prevedere i valori futuri.
Qualsiasi serie temporale "non stagionale" che mostri schemi e non sia un white noise può essere modellata con i modelli ARIMA.
Un modello ARIMA è caratterizzato da 3 termini: p, d, q
p è l'ordine del termine AR
q è l'ordine del termine MA
d è il numero di differenze necessarie per rendere stazionaria la serie temporale
Se una serie temporale presenta modelli stagionali, è necessario aggiungere termini stagionali e diventa SARIMA, abbreviazione di "ARIMA stagionale".
Il primo passo per costruire un modello ARIMA è rendere stazionarie le serie temporali; questo perché il termine "Regressione automatica" in ARIMA significa che si tratta di un modello di regressione lineare che utilizza i propri ritardi come predittori. I modelli di regressione lineare funzionano meglio quando i predittori non sono correlati e sono indipendenti l'uno dall'altro.
Per rendere stazionaria la serie, l'approccio più comune è quello di differenziarla. Cioè, sottrai il valore precedente dal valore corrente. A volte, a seconda della complessità della serie, potrebbe essere necessaria più di una differenziazione.
Il valore di d, quindi, è il numero minimo di differenze necessarie per rendere stazionaria la serie. E se la serie temporale è già stazionaria, allora d = 0.
Per quanto riguarda gli altri termini, "p" è l'ordine del termine "Auto Regressive" (AR). Si riferisce al numero di ritardi di Y da utilizzare come predittori; "q" è l'ordine del termine "Media mobile" (MA). Si riferisce al numero di errori di previsione ritardati che dovrebbero entrare nel modello ARIMA.
TROVARE L'ORDINE d¶
Lo scopo di differenziare la serie è quello di renderla stazionaria.
Differenziare troppo la serie potrebbe essere una cosa negativa, in quanto una serie troppo differenziata può essere nuovamente stazionaria.
Il giusto ordine di differenziazione è la differenza minima richiesta per ottenere una serie quasi stazionaria che si aggira intorno a una media definita e il grafico ACF raggiunge lo zero abbastanza rapidamente.
Se le autocorrelazioni sono positive per molti numeri di ritardi (10 o più), la serie necessita di ulteriori differenze. D'altra parte, se la stessa autocorrelazione del ritardo 1 è troppo negativa, la serie è probabilmente eccessivamente differenziata.
Per prima cosa si verifica che la serie è stazionaria, usando l'Augmented Dickey Fuller Test. Nel caso in cui la serie sia stazionaria, non sarà necessaria la differenziazione e d sarà uguale a 0.
L'ipotesi nulla del test ADF è che la serie temporale non sia stazionaria. Quindi, se il valore p del test è inferiore al livello di significatività (0,05), si rifiuta l'ipotesi nulla e si deduce che la serie temporale è effettivamente stazionaria.
Per calcolare l'ACF elimino i valori nulli, considerando solo i valori numerici.
from statsmodels.tsa.stattools import adfuller
from numpy import log
# ts_na = df_grouped_date[['Date', 'Minuti ciclo']]
result = adfuller(ts_f['Minuti ciclo'].dropna())
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
Il p-value è maggiore del livello di significatività . A questo punto, procedo con la differenziazione.
import matplotlib.pyplot as plt
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
# pyplot.figure(figsize=(12,5), dpi=100)
# plot_acf(ts_na['Minuti ciclo'].dropna(), lags = 200)
# pyplot.show()
# plot_acf(ts_na['Minuti ciclo'].diff().dropna(), lags = 200)
# pyplot.show()
# plot_acf(ts_na['Minuti ciclo'].diff().diff().dropna(), lags = 200)
#Â pyplot.show()
#print if you want to visualize
#print(dta.head())
fig, ax = plt.subplots(1,3,figsize=(20,5))
sm.graphics.tsa.plot_acf(ts_f['Minuti ciclo'].dropna(), lags=60, ax=ax[0])
sm.graphics.tsa.plot_acf(ts_f['Minuti ciclo'].diff().dropna(), lags=60, ax=ax[1])
sm.graphics.tsa.plot_acf(ts_f['Minuti ciclo'].diff().diff().dropna(), lags=60, ax=ax[2])
plt.show()
Raggiunge la stazionarietà con 1 differenza. Infatti con 2 differenze i primi lag sono molto più marcati, quindi considero solo 1 differenza!
result = adfuller(ts_f['Minuti ciclo'].diff().dropna())
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
Ora è differenziata
TROVARE L'ORDINE p¶
Il passaggio successivo consiste nell'identificare se il modello necessita di termini AR. Per scoprirlo, bisogna controllare il grafico PACF.
L'autocorrelazione parziale può essere immaginata come la correlazione tra la serie e il suo ritardo, dopo aver escluso i contributi dai ritardi intermedi. Quindi, PACF trasmette in qualche modo la pura correlazione tra un ritardo e la serie.
Qualsiasi autocorrelazione in una serie stazionaria può essere rettificata aggiungendo un numero sufficiente di termini AR. Quindi, inizialmente assumiamo che l'ordine del termine AR sia uguale al numero di ritardi che supera il limite di significatività nel grafico PACF.
import matplotlib.pyplot as pyplot
plot_pacf(ts_f['Minuti ciclo'].diff().dropna(), lags = 60)
pyplot.show()
il primo lag risulta significativo, come il secondo. Il terzo non risulta molto significativo, quindi fisso il valore di p a 2
TROVARE L'ORDINE q¶
Come per l'ordine p bisogna guardare la PACF, per definire l'ordine q (MA) bisogna visualizzare la ACF.
plot_acf(ts_f['Minuti ciclo'].diff().dropna(), lags = 60)
pyplot.show()
Il primo ed il secondo lag sono fuori dagli intervalli di confidenza, quindi considero q = 2.
MODELLO ARIMA¶
Ora che sono stati definiti i valori p, d, q, viene implementato il modello ARIMA.
from statsmodels.tsa.arima_model import ARIMA
# 2,1,2 ARIMA Model
model1 = ARIMA(ts_f['Minuti ciclo'], order=(6,1,2))
model1_fit = model1.fit(disp=0)
print(model1_fit.summary())
import matplotlib.pyplot as plt
# Plot residual errors
residuals = pd.DataFrame(model1_fit.resid)
fig, ax = plt.subplots(1,2)
residuals.plot(title="Residuals", ax=ax[0])
residuals.plot(kind='kde', title='Density', ax=ax[1])
plt.show()
# plot real vs fitted
model1_fit.plot_predict(dynamic=False)
plt.show()
fc_train_series1 = pd.Series(model1_fit.predict(typ='levels'), index=train.index)
fc_train_series1 = fc_train_series1.fillna(fc_train_series1.mean())
# RMSE ON TRAINING
from sklearn.metrics import mean_squared_error
from math import sqrt
rmse1_train = sqrt(mean_squared_error(train, fc_train_series1))
rmse1_train
La tabella dei coefficienti ci dice che il coefficiente dell'AR6 è poco significativo ed il p value (colonna P>|z|) è al di fuori dell'intervallo di confidenza per MA1 ed MA2. Provo un ARIMA(6,1,3)
# 6,1,3 ARIMA Model
model2 = ARIMA(ts_f['Minuti ciclo'], order=(6,1,3))
model2_fit = model2.fit(disp=0)
print(model2_fit.summary())
L'AIC del modello è leggermente diminuito ed è una cosa positiva. I coefficienti di AR ed MA sono significativi per AR1 e MA2, mentre MA1 > 0.05 e quindi non è significativo.
Procedo con la stampa dei residui
import matplotlib.pyplot as plt
# Plot residual errors
residuals = pd.DataFrame(model2_fit.resid)
fig, ax = plt.subplots(1,2)
residuals.plot(title="Residuals", ax=ax[0])
residuals.plot(kind='kde', title='Density', ax=ax[1])
plt.show()
La media è vicina allo zero e la varianza è uniforme.
# plot real vs fitted
model2_fit.plot_predict(dynamic=False)
plt.show()
fc_train_series2 = pd.Series(model2_fit.predict(typ='levels'), index=train.index)
fc_train_series2 = fc_train_series2.fillna(fc_train_series2.mean())
# RMSE ON TRAINING
from sklearn.metrics import mean_squared_error
from math import sqrt
rmse2_train = sqrt(mean_squared_error(train, fc_train_series2))
rmse2_train
Una volta identificato il modello, ora effettuiamo la previsione dei valori futuri.
FORECASTING¶
from statsmodels.tsa.stattools import acf
# Create Training and Test
# train = ts_f['Minuti ciclo'][:366]
# test = ts_f['Minuti ciclo'][367:571]
# Build Model
model = ARIMA(train, order=(1, 1, 2))
fitted = model.fit(disp=-1)
# Forecast
fc, se, conf = fitted.forecast(len(test), alpha=0.05) # 95% conf
# Make as pandas series
fc_series = pd.Series(fc, index=test.index)
lower_series = pd.Series(conf[:, 0], index=test.index)
upper_series = pd.Series(conf[:, 1], index=test.index)
# Plot
plt.figure(figsize=(12,5), dpi=100)
plt.plot(train, label='training')
plt.plot(test, label='actual')
plt.plot(fc_series, label='forecast')
plt.fill_between(lower_series.index, lower_series, upper_series,
color='k', alpha=.15)
plt.title('Forecast vs Actuals')
plt.legend(loc='upper left', fontsize=8)
plt.show()
from sklearn.metrics import mean_squared_error
from math import sqrt
rms1 = sqrt(mean_squared_error(test, fc_series))
rms1
La previsione senza la componente stagionale non riesce ad identificare correttamente le previsioni.
Introduco la componente stagionale.
SARIMA¶
Inizialmente provo l'auto.arima, poi sviluppo il modello.
Provo l'auto.arima
import pmdarima as pm
# Seasonal - fit stepwise auto-ARIMA
smodel = pm.auto_arima(train, start_p=1, start_q=1,
test='adf',
max_p=6, max_q=3, m=6,
start_P=0, seasonal=True,
d=None, D=1, trace=True,
error_action='ignore',
suppress_warnings=True,
stepwise=True)
smodel.summary()
l'auto arima identifica la previsione migliore con un ARIMA(1,0,1)(1,1,1) di periodo 6
# sARIMA Model
sarima_model = SARIMAX(ts_f['Minuti ciclo'], order=(1,0,1), seasonal_order = (1,1,1,6))
sarima_fit = sarima_model.fit(disp=0)
# plot real vs fitted
sarima_fit.plot_diagnostics()
plt.show()
fc_train_series_smodel = pd.Series(sarima_fit.predict(typ='levels'), index=train.index)
fc_train_series_smodel = fc_train_series_smodel.fillna(fc_train_series_smodel.mean())
# RMSE ON TRAINING
from sklearn.metrics import mean_squared_error
from math import sqrt
rmse_smodel_train = sqrt(mean_squared_error(train, fc_train_series_smodel))
rmse_smodel_train
# Forecast
n_periods = len(test)
fitted, confint = smodel.predict(n_periods=n_periods, return_conf_int=True)
index_of_fc = pd.date_range(ts['Minuti ciclo'].index[-1], periods = n_periods, freq='D')
# make series for plotting purpose
fitted_series = pd.Series(fitted, index=test.index)
lower_series = pd.Series(confint[:, 0], index=test.index)
upper_series = pd.Series(confint[:, 1], index=test.index)
# Plot
plt.figure(figsize=(12,5), dpi=100)
plt.plot(train, label='training')
plt.plot(test, label='actual')
plt.plot(fitted_series, label='forecast', color='darkgreen')
plt.fill_between(lower_series.index,
lower_series,
upper_series,
color='k', alpha=.15)
plt.fill_between(lower_series.index, lower_series, upper_series,
color='k', alpha=.15)
plt.title("SARIMA - Final Forecast")
plt.legend(loc='upper left', fontsize=8)
plt.show()
# calcolo residui
residuals = fitted_series - test['Minuti ciclo']
plt.figure(figsize=(10,4))
plt.plot(residuals)
plt.axhline(0, linestyle='--', color='k')
plt.title('Residuals from SARIMA Model', fontsize=20)
plt.ylabel('Error', fontsize=16)
Le previsioni sono abbastanza buone, si intravede un trend crescente e viene identificata correttamente la stagionalità con picchi negativi al sabato.
from sklearn.metrics import mean_squared_error
from math import sqrt
rmse_sarima = sqrt(mean_squared_error(test, fitted_series))
rmse_sarima
ts_f
# save
ts_f.to_csv(folder_risultati + '//df_min_ciclo_nosunday_no0.csv', decimal = ',', index = True)
SARIMA WITH FOURIER TERMS¶
# r script process data and generate results
# Load fourier.csv, original file
pre_fourier = pd.read_csv(folder_risultati+"/prediction_fourier.csv", sep = ",").replace('.', ',', regex=True).astype(float)
lower_fourier = pd.read_csv(folder_risultati+"/lower_fourier.csv", sep = ",").replace('.', ',', regex=True).astype(float)
upper_fourier = pd.read_csv(folder_risultati+"/upper_fourier.csv", sep = ",").replace('.', ',', regex=True).astype(float)
# Forecast
test_2 = test[:-2]
n_periods = len(test_2)
index_of_fc = pd.date_range(ts['Minuti ciclo'].index[-1], periods = n_periods, freq='D')
# make series for plotting purpose
fitted_series = pd.Series(pre_fourier['x'].to_numpy(), index=test_2.index)
lower_series = pd.Series(lower_fourier['x'].to_numpy(), index=test_2.index)
upper_series = pd.Series(upper_fourier['x'].to_numpy(), index=test_2.index)
# Plot
plt.figure(figsize=(12,5), dpi=100)
plt.plot(train, label='training')
plt.plot(test, label='actual')
plt.plot(fitted_series, label='forecast', color='darkgreen')
plt.fill_between(lower_series.index,
lower_series,
upper_series,
color='k', alpha=.15)
plt.fill_between(lower_series.index, lower_series, upper_series,
color='k', alpha=.15)
plt.title("ARIMA - Forecast with fourier terms")
plt.legend(loc='upper left', fontsize=8)
plt.show()
rmse_arima_fourier = sqrt(mean_squared_error(test_2, fitted_series))
rmse_arima_fourier